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FOREWORD 


The  performance  testing  of  the  NAVPAC  multisatellite  receiver  was  accomplished  by  observing 
passes  of  particular  satellites  of  the  Navy  Navigation  Satellite  System.  These  passes  were  stored  on  magnetic 
tape  and  sent  to  the  Naval  Surface  Weapons  Center,  Dahlgren  Laboratory  (NSWC/DL)  for  reduction.  This 
report  describes  the  formulation  which  comprised  this  data  processor,  and  which  was  current  in  June  1976. 

This  formulation  is  a small  part  of  the  entire  NAVPAC  project  at  NSWC/DL.  The  instrument  itself 
was  designed  and  built  at  the  Applied  Physics  Laboratory  under  the  direction  of  R.  E.  Willison.  At 
Dahlgren  the  technical  direction  and  project  management  were  by  Mr.  E.  Swift  with  programming  support 
by  Mr.  A.  Fisher. 


Released  by; 


R.  A.  Niemann 

Head,  Warfare  Analysis  Department 
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INTRODUCTION 


The  NAVPAC  multisatellite  receiver  has  the  capability  of  receiving  simultaneously  three  dual- 
frequency  navigation  signals  transmitted  from  the  Navy  Navigation  Satellite  System.  The  NAVPAC  receiver 
was  developed  by  the  Applied  Physics  Laboratory  of  the  Johns  Hopkins  University  and  consists  of  the  fol- 
lowing subsystems;  1 50/400-MHz  antenna,  the  multisatellite  receiver,  a stable  reference  oscillator  and 
time  code  generator,  a data  processor  and  core  memory  for  storage  of  the  data,  and  a tape  recorder  con- 
troller for  final  disposition  of  the  accumul.ited  data.  This  report  describes  the  formulation  used  to  convert 
the  observations  of  Doppler  shifted  NAVSAT  transmissions  during  one  satellite  pass  into  a two-dimensional 
navigation  solution  for  the  position  of  the  NAVPAC  antenna. 

When  the  NAVPAC  tape  containing  Doppler  observations  arrives  at  the  processing  site,  the  first  step 
of  reduction  (Figure  1)  involves  unpacking  the  data  words  and  separating  them  into  their  constituent  parts. 
Table  1 illustrates  the  composition  of  each  of  the  four  types  of  words.  The  Doppler  time  mark  word  (DTM) 
is  composed  of  the  receiver  identification  (RID)  and  the  NAVPAC  time.  The  Doppler  word  (DW)  contains 
a receiver  identification,  a raw  Doppler  count,  a refraction  count,  and  the  corresponding  time  over  count. 
The  Doppler  and  refraction  counts  will  be  described  in  detail  in  a later  section.  The  two-minute  time  mark 
word  (TMTM)  has  the  NAVSAT  identification  word  (NID),  the  receiver  identification,  and  the  NAVPAC 
time  at  the  instant  of  reception  of  the  two-minute  mark.  The  vehicle  time  mark  (VTM)  is  the  final  word 
type.  Any  subset  of  the  above  information  can  be  printed  as  required  by  the  various  input  options. 

The  second  step  in  the  processing  is  the  correction  of  the  NAVPAC  clock  by  comparison  to  the  time 
signals  transmitted  by  the  navigation  satellites  and  recorded  as  NAVPAC  times  in  the  two-minute  time 
mark  words.  The  corrections  obtained  are  added  to  the  Doppler  time  marks  in  order  to  establish  the  pre- 
cise (<50-microsecond  error)  time  reference  for  a particular  Doppler  count. 

The  third  operation  is  to  form  passes  from  the  chronological  collection  of  data.  A satellite  pass  is 
defined  to  be  the  duration  of  time  during  which  a particular  satellite  is  “visible”  from  the  receiving  station. 
For  this  operation  a trajectory  tape,  containing  predicted  (or  best  fit  if  reduction  is  done  at  a later  time) 
earth  fixed  positions  of  all  satellites  of  interest,  is  needed  as  a key  to  establish  the  beginnings  and  ends  of 
passes.  The  trajectory  is  especially  needed  to  mark  the  ends  of  passes  since  NAVPAC  gives  no  definitive 
indication  when  it  loses  lock  on  a satellite.  A related  function  is  to  identify  short,  low  elevation,  passes 
which  may  not  contain  a TMTM  between  rise  and  set.  Such  passes  cannot  be  identified  from  NAVPAC 
data  alone  since  the  satellite  identification  is  present  only  in  TMTM’s  (Table  I).  Finally,  the  trajectory  tape 
is  used  to  establish  which  data  will  be  processed.  If  a particular  satellite  is  not  on  the  trajectory  file,  it  will 
not  be  considered  further. 


Table  1.  Composition  of  NAVPAC  Doppler  Data  Words 


Type  of  Word 

RID 

NID 

Doppler 

Count 

Time  Over 
Count 

Refraction 

Count 

NAVPAC 

Clock 

DTM 

X 

X 

DW 

X 

X 

X 

X 

TMTM 

X 

X 

X 

VTM 

X 

X 

1 


The  data  which  is  accepted  from  each  pass  is  written  in  time  order  in  a format  making  it  compatible 
with  the  CELEST  orbit  computation  program.*  Alternately,  if  the  full  capability  of  CELEST  is  not 
required,  a smaller  processor  can  be  used  as  the  final  step.  This  program,  which  uses  the  data  to  produce  a 
two-dimensional  navigation  on  a pass-by-pass  basis,  is  described  in  later  sections. 


PASS  RECOGNITION 

All  Navy  navigation  satellites  are  synchronized  to  Coordinated  Universal  Time  (UTC)  to  within  about 
50  microseconds.  Therefore  the  two-minute  time  marks  received  at  NAVPAC  provide  a UTC  time  reference 
having  a resolution  of  about  this  accuracy.  Use  is  made  of  these  time  marks  to  establish  the  beginning  and 
end  time  of  each  Doppler  count.  Each  pass  begins  witli  a Doppler  time  mark,  indicating  that  a receiver  has 
successfully  acquired  a satellite  signal.  The  time  indicated  by  this  Doppler  time  mark  is,  however,  not  the 
time  at  which  the  count  began,  but  rather  the  time  30  NAVPAC  seconds  later.  This  procedure  was  adopted 
to  eliminate  a multitude  of  Doppler  time  marks  which  would  have  occurred  at  every  false  acquisition  signal. 
The  30-second  delay  requires  that  NAVPAC  remain  locked  onto  a satellite  signal  for  30  seconds  before  such 
a fact  is  acknowledged  by  the  data  processor. 

The  first  Doppler  word  following  the  Doppler  time  mark  contains  the  Doppler  count  for  the  period 
beginning  at  signal  acquisition  and  ending  at  the  time  indicated  by  the  sum  of  the  time  over  count  and  the 
time  in  the  Doppler  time  mark  word.  Subsequent  Doppler  words  are  produced  at  each  succeeding  observa- 
tion. Observations  occur  at  intervals  of  30  NAVPAC  seconds  plus  time  over  counts.  A NAVPAC  unit  of 
time  is  approximately  equal  to  1 .2  X 10~^  UTC  seconds;  more  precisely,  2.5  X 10^  NAVPAC  time  units 
equal  exactly  30  NAVPAC  seconds.  Time  over  counts  are  always  less  than  128  microseconds.  After  eiglit 
consecutive  Doppler  words,  another  Doppler  time  mark  is  produced.  The  time  between  successive  Doppler 
time  marks  should  equal  240  seconds  plus  the  sum  of  the  intervening  time  over  counts.  Interspersed  with 
the  Doppler  time  marks  and  Doppler  words  are  the  two-minute  time  marks  indicating  the  NAVPAC  time 
at  each  even  two-minute  interval. 

An  example  of  a four-minute  segment  from  one  satellite  is  shown  in  Table  2.  The  time  between  the 
two  Doppler  time  marks  is  240.000162  seconds,  and  the  sum  of  the  time  over  counts  for  the  intervening 
eight  Doppler  words  is  162  microseconds,  indicating  agreement.  Figure  2 is  the  same  example  in  a different 
format.  Figure  2 illustrates  that  the  DTM’s  occur  immediately  before  the  ith  and  (i  + 8)th  Doppler  words 
(where  i = 1 in  Figure  2),  and  includes  the  eight  time  over  counts.  Note  that  the  beginning  of  the  pass  is 
actually  at  3467.274860  seconds  and  therefore  the  first  TMTM  occurred  after  satellite  acquisition. 

Since  NAVPAC  has  the  capability  of  receiving  signals  from  three  satellites  simultaneously,  the  data 
from  each  satellite  can  be  intermixed  with  the  other  two.  Also,  the  possibility  exists  that  receivers  and 
satellites  may  switch  without  warning.  For  these  reasons  the  pass  recognition  and  sorting  logic  become 
very  complex.  The  Doppler  words  give  no  direct  information  as  to  which  satellite  the  present  count  refers, 
only  the  receiver  identification.  The  receiver  identification  in  turn  is  related  to  a satellite  through  informa- 
tion supplied  by  the  TMTM’s.  The  TMTM  decoder  cycles  through  the  active  receivers  in  sequence.  Thus  if 
all  are  tracking,  a gap  of  six  minutes  can  occur  before  it  is  confirmed  that  each  receiver  is  still  tracking  the 
same  satellite.  If  a receiver  has  switched,  then  the  DTM’s  may  indicate  a time  difference  that  is  not  240  + 
time  over  seconds  somewhere  within  the  previous  sbe  minutes.  Should  this  be  the  case,  the  number  of 
possibilities  of  what  anomaly  may  have  occurred  is  quite  large  and  unfortunately  not  unique.  Some  of  the 
anomalies  that  can  occur  and  which  arc  recognized  are  the  following: 
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Seconds 


Figure  2.  Graphical  Representation  of  Table  2 
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Table  2.  Representative  Doppler  Data  Beginning  a New  Pass 


Record  Type 

Time 

Time  Over  Count 

Raw  Doppler 
Count 

NAVPAC 

Units 

Ai  seconds 

NAVPAC  Units 

Seconds 

TMTM 

2900728587 

3481.114304 

DTM 

2914395717 

3497.274860 

DWi 

30 

36.0 

710435 

DWj 

1 

1.2 

710848 

DW3 

10 

12.0 

711591 

DW4 

22 

26.4 

712731 

TMTM 

3000926143 

3601.111372 

DW5 

7 

8.4 

714363 

DWfi 

19 

22.8 

716616 

DW7 

25 

30.0 

719661 

DWg 

21 

25.2 

723735 

TMTM 

3100923875 

3721.108650 

DTM 

3114395852 

3737.275022 

a.  Loss  of  Lock.  A receiver  once  locked  and  tracking  a NAVSAT  may  lose  lock  at  some  point  in  the 
pass.  It  is  possible  that  this  signal  may  be  reacquired  by  another  receiver  (which  was  in  the  search  mode) 
before  it  is  found  again  by  the  original  tracking  receiver.  In  this  case  the  characteristic  symptom  is  that  the 
NID  will  switch  to  a new  RID. 

b.  Lock  on  Bogie.  A receiver  once  locked  and  tracking  a NAVSAT  may  be  enticed  to  capture  a 
strong  but  extraneous  signal  not  having  the  required  signal  structure.  The  NAVPAC  tracking  logic  will 
notice  both  the  lack  of  correct  structure  and  presumably  little  or  no  Dopple; . The  signal  will  be  dropped, 
but  in  the  meantime  it  has  broken  the  continuity  of  the  pass. 

c.  High-Priority  Switching.  Two  receivers,  I and  II,  are  initially  locked  onto  two  satellite  signals:  a 
strong  signal  “s”  and  a weak  signal  “w.”  Receiver  I is  tracking  s and  II  is  tracking  w.  Suppose  the  NAVPAC 
priority  logic  has  given  II  priority  over  I.  The  priority  logic  decides  that,  if  two  receivers  acquire  and  lock 
onto  the  same  satellite,  the  one  that  acquired  first  will  continue,  but  the  later  receiver  will  be  unlocked  and 
resume  its  search  mode.  Suppose  the  two  signals  are  such  that  the  Doppler  frequencies  coincide  at  some 
point  in  the  pass.  It  is  possible  that  the  high-priority  receiver  (II)  may  jump  to  the  stronger  signal  s without 
losing  lock  or  otherwise  indicating  what  happened.  Since  II  still  has  priority,  the  NAVPAC  logic  will  kick  I 
off  s when  it  discovers  that  two  receivers  are  tracking  the  same  satellite.  Thus  the  correct  pass  structure  for 
s begins  on  receiver  1 and  ends  on  II,  with  the  pass  for  w on  II  being  incomplete. 

d.  Satellite  Swapping.  This  case  is  similar  to  the  above  except  that  now  receiver  1 acquires  the  weak 
signal  within  the  same  30-second  interval  that  II  switches  to  the  strong  signal.  For  this  to  happen,  the  signal 
strengths  must  be  nearly  equal  and  the  Doppler  slopes  nearly  identical. 
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NAVPAC  TIME  CORRECTION 


To  establish  the  correspondence  between  UTC,  designated  by  the  symbol  t (seconds),  and  NAVPAC 
time  T (in  NAVPAC  units),  the  reset  time  t,  must  be  correct  to  within  one  minute.  The  reset  time  is  defined 
as  the  UTC  instant  when  NAVPAC  clock,  via  a command,  is  reset  to  zero.  The  first  TMTM  received  after 
resetting  can  then  be  used  to  refine  the  estimate  of  tr  to  an  accuracy  better  than  one  second.  The  procedure 
is  as  follows:  The  time  (seconds)  between  reception  of  the  first  TMTM  and  the  reset  time  is  1.2  X I0~®  t. 
The  UTC  of  the  TMTM  is  therefore  tt  = t,  + 1 .2  X 10“®  t ± one  minute,  but  it  is  known  that  all  TMTM’s  are 
transmitted  on  the  even  minute.  Therefore  a better  approximation  to  the  reset  time  is  obtained  by  round- 
ing t|  to  the  nearest  even  minute  tj  and  calculating  a new  reset  time  t|^ 

tR  = tj  - 1.2  X 10-6  tj  = tjlfounded  • 

It  is  inevitable  that  NAVPAC  time  will  wander  slightly  from  the  UTC  standard.  In  order  to  correct 
for  this,  an  expression  relating  the  two  is  needed.  The  expression  should  be  such  that  given  any  t,  the 
corresponding  t can  be  evaluated  to  within  ± 50 /iseconds.  The  relation  chosen  is  a polynomial: 

t(T)  = tg  + a + (1  + b)T  + CT^  (I) 

The  form  of  this  polynomial  models  the  initial  offset  by  t(^  + a,  a rate  error  by  I + b,  and  the  aging 
phenomena  by  c.  The  three  parameters  a,  b,  c are  small  corrections  to  the  values  t|^,  one,  and  zero  expected 
from  a clock  which  follows  UTC  exactly. 

All  TMTM’s  in  a preset  interval  are  used  to  evaluate  the  parameters  a,  b,  c in  a least-squares  solution. 
The  time  of  reception  of  the  ith  TMTM  is  the  observed  quantity  T^.  The  UTC  corresponding  to  this  instant 
is 

‘i  = tMi  + tpi  + td.  (2) 

where  t|^j  is  the  time  of  transmission  of  the  TMTM  signal  at  the  satellite  (an  even  minute),  tpj  is  the  propaga- 
tion time  from  satellite  to  NAVPAC,  and  t^  is  a decode  delay  in  NAVPAC.  The  value  for  tpj  is  obtained 
from  the  trajectory  file  (Figure  1)  by  calculating  the  slant  range  at  each  tj.  Thus  tpj  = l/c  (slant  range  at  tj). 
The  error  introduced  in  tpj  by  not  calculating  the  slant  range  correctly  and  by  not  including  refraction  is 
less  than  a microsecond.  This  error  will  be  ignored.  The  difference  between  Equations  (1)  and  (2)  is  the 
error  which  is  minimized  in  the  least-squares  sense: 

Cj  = a + (1  + b)Tj  + CT?  - tMj  - tpj  - tj  + t,^  (3) 

The  decode  delay  t^  is  inserted  in  Equation  (2)  because  this  delay  affects  only  TMTM’s,  and  its  presence 
must  be  eliminated  in  order  to  obtain  a correct  value  for  a in  Equation  (1). 

The  consistency  of  the  solution  is  checked  by  inserting  the  observations  Tj  into  Equation  ( 1)  and 
calculating  the  corresponding  tj’s.  If  these  calculated  values  differ  from  the  corrected  TMTM’s  obtained 
from  Equation  (2)  by  a preset  tolerance,  the  observation  is  rejected.  A new  solution  is  then  performed 
using  only  the  observations  which  pass  the  above  test.  This  iteration  process  eliminates  from  consideration 
TMTM’s  which  do  not  agree  with  the  majority.  Equation  (1)  can  now  be  used  to  convert  any  NAVPAC 
time  T to  its  UTC  equivalent. 
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SYNTHESIS  OF  THE  RANGE  DIFFERENCE  OBSERVATION 

The  NAVPAC  receivers  are  designed  to  acquire  and  track  dual-frequency  Doppler  signals  from  the 
NAVSATS  and  convert  them  to  Doppler  counts  and  refraction  counts  at  intervals  of  about  30  seconds. 

Since  the  relative  velocity  between  NAVPAC  and  NAVSAT  is  continuously  varying,  the  period  of  the  Dop- 
pler cycles  continuously  varies.  As  illustrated  in  Figure  3,  the  time  over  count  enters  from  the  desire  to 
keep  the  number  of  cycles  (counts)  in  an  observation  an  integer.  This  is  accomplished  by  recording  the 
time  in  excess  of  30  seconds  needed  to  complete  the  cycle  which  was  in  progress  at  the  30-second  mark. 
Thus  an  integral  number  of  Doppler  cycles,  the  Doppler  count,  occurs  in  the  interval  (30  + time  over) 
seconds. 

The  signal  processing  which  NAVPAC  performs  changes  the  two  received  satellite  frequencies  and 
1^2  into  a Doppler  count  n^  and  a refraction  count  nj.  In  order  to  explain  the  process,  the  following  parame- 
ters need  to  be  defined: 

f|  = 150  X 10^  Hz,  the  nominal  frequency  of  the  transmitted  signal 

fg  = - 1 2000  Hz,  the  frequency  offset 

fb  = an  unknown  time  variable  bias,  NAVSAT  dependent 

fp  = an  unknown  time  variable  bias  in  NAVPAC 

q =8/3 

fbi  = the  ionospheric  plasma  frequency 

For  solution  purposes,  the  two  biases  fj,  and  fp  are  assumed  to  be  constant  over  each  20-minute 
interval  comprising  a pass.  Each  NAVSAT  transmits  two  frequencies:  (fj  + fq  + fb)  and  q(f j + f q + fb). 

The  frequencies  received  at  NAVPAC  which  correspond  to  these  have  added  to  them  a Doppler  component 
and  an  ionospheric  component.  These  are  modeled  in  Equation  (4): 

>^l=(f,  + fo  + fb)  + (4a) 

''2  = q(fl  +fo  + fb)  + (4b) 

The  sign  convention  used  for  the  relative  velocity  v is  that  v < 0 if  the  transmission  path  length  is  decreasing. 
In  Equations  (4),  the  earth’s  magnetic  field  is  not  considered  and  so  the  index  of  refraction  n,  and  the 
received  frequency  u are2.3 

«'  = (fi  + fq  fb)  • 

If  n is  expanded  as  a binomial  series,  the  result  is 

n - 1 -y  f2(f,  + t'o  ^ f,)-2  -±f4(f,  + + f^)-4  + ... 

The  high-order  terms  are  normally  negligible  and  so  are  ignored  in  Equations  (4). 
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Figure  3.  Diagram  of  the  Elements  of  a Doppler  Count 


In  NAVPAC,  the  Doppler  frequency  is  obtained  by  subtracting  P2  from  a locally  generated  frequency 
q(fl  + fo  + fp).  In  symbolic  form,  the  result  is 


+fo  + <‘n)-''2=q(fn-fb)  + ‘l— (f]  +''o  + fb) 


The  Doppler  count  is  obtained  by  integrating  this  Doppler  frequency  over  the  observation  interval  (tj  - tj.) ). 
This  is  equivalent  to  counting  the  number  of  Doppler  cycles. 


nd=  f lq(fi  + fo  + fn)-‘'2l 

■''i-l 

-iJ 


dt 


"d  = q(fn  - fb) 


• f2  q 


rl 


(fi+fo  + fb)-^fi5  V(fi+fo  + fb)' 


-11: 


vdt 


(S) 


The  refraction  frequency  is  obtained  by  subtracting  t'j  from  q~^«'2- 

q'‘*'2  - "1  ^0  + fb)'^(q'^  - 1) 

Integrating  this  over  the  same  time  interval  results  in  the  refraction  count  n^, 


nr  = 4^Y)f^-^^^l  +fo  + fbr'(q‘^-l)£  vdt 


(6) 


In  Equation  (6),  the  factor  4 is  inserted  by  NAVPAC.  If  n,  is  multiplied  by  an  appropriate  constant  the 
result  can  be  used  to  simplify  Equation  (5).  This  form  is  written  as  Equation  (7): 


q~'  Vr  1 

4 


-1 


.q'^  - 1 


+ <'0  + %)-' /*' 
‘'•i-l 


vdt 


(7) 


The  constant  was  chosen  so  that  the  last  term  on  the  right  in  Equation  (5)  is  the  same  as  in  Equation  (7). 

If  these  two  expressions  arc  summed,  the  effects  of  refraction  are  eliminated  from  the  Doppler  count.  The 
corrected  Doppler  count  is  therefore  n^: 


q”*  6 


(8) 


These  raw  counts  are  now  converted  to  range  differences  so  they  can  be  used  as  observations  in  the 
remainder  of  the  formulation.  Performing  the  integration  in  Equation  (5),  where  p is  the  range  from  NAV- 
SAT  to  NAVPAC,  gives 


"c  = q(fn  - fbXtj  - tj_, ) + Y (f|  + fo  )1  > 


9 


where 


dt 


Finally,  the  observation  expression  is 


O: 


cn. 


q(fl  +l'o  + 'b> 


c(fn-fb) 

-—7-77-  (ti-ti_,)  + p(tj)-p(ti.,) 

1 ^ In  ^ ‘K 


(9) 


which  is  in  units  of  meters  and  is  the  range  difference  expression  Oj,  for  the  NAVPAC  corrected  count  n^ 
ending  at  time  tj.  Wlien  Fquation  (q)  is  evaluated,  fn  is  taken  to  be  zero. 


MODEL  FORMULATION 


The  primary  result  desired  from  this  formulation  is  an  indication  of  the  level  of  performance  possible 
with  the  NAVPAC  hardware  compared  with  the  performance  of  other  comparable  systems.  To  this  end, 
this  formulation  must  be  similar  to  that  used  to  process  satellite  data  on  a routine  basis.  Therefore  many  of 
the  techniques  described  here  were  chosen  because  of  their  use  in  the  CELEST  orbit  computation  program.^ 

The  NAVPAC  antenna  location  can  be  expressed  in  either  earth-fixed  Xj,,  yg,  Zg  coordinates  or  in 
geodetic  latitude  0g,  longitude  Xg,  and  height  hg.  Given  an  ellipsoid  of  revolution  about  the  z axis  with 
eccentricity  e and  semimajor  axis  a,  the  relationship  between  the  two  sets  of  coordinates  is 


(0  + hg)cos0g  cosXg 

(10a) 

(Q  + hg)  cos0g  sin  Xg 

(10b) 

(Q  + hg ) sin  0g  - Qe^  sin  i/ig 

(10c) 

where 


0 = a|  I - sin^  0g] 

In  the  formulation  that  follows,  earth  fixed  coordinates  are  used;  however,  in  many  cases  the  station 
location  is  given  as  geodetic  coordinates  on  an  ellipsoid.  Equations  (10)  are  included  here  to  indicate  how 
the  transformation  is  accomplished. 

The  previous  section  described  how  range  difference  information  is  obtained  from  the  Doppler 
observations.  The  value  of  range  difference  at  each  observation  time  tj  is  Oj.  The  corresponding  value 
calculated  from  information  residing  on  the  trajectory  file  is  Cj.  The  difference  of  these  two  (Oj  - Cj)  is 
the  error  which  is  to  be  minimb.ed  iti  a least-squares  sense. 

To  form  the  expression  for  the  computed  range  difference,  the  range  from  satellite  to  NAVPAC  is 
necessary  at  the  instant  of  transmission  at  the  beginning  and  end  of  a Doppler  count  interval.  Knowledge  of 
the  propagation  time  is  needed  to  obtain  both  values.  The  trajectory  file  provides  estimates  of  the 
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NAVSAT  earth  fixed  position  x,,  y,,  Zj  at  equal  intervals  in  time,  but  positions  at  other  times  are  available 
by  use  of  eight-point  LaGrangian  interpolation. 

In  order  to  establish  the  propagation  time  at  the  beginning  and  end  of  a Doppler  count,  an  iterative 
prtKedure  is  used  as  follows.  The  slant  range  from  NAVPAC  to  NAVSAT  is  the  magnitude  of  the  difference 
of  their  respective  vectors  7^  and  Fj.  In  earth  fixed  coordinates,  the  vector  to  NAVPAC  is  Fg,  and  to  NAV- 
SAT u- 


h = XsX  + yS  + ZjZ 

Fg  = xgx  + ygy  + zgz 


(lla) 


The  range  is  the  magnitude  of  p 

P = ts  - fp  = (x,  - Xg)x  + (ys  - yg)y  + (z^  - Zg)z  ( 1 2) 

The  proper  propagation  time  is  obtained  when  Fj  is  evaluated  at  the  instant  of  transmission  of  the  Doppler 
signal  tj^^j,  while  Fg  is  evaluated  at  the  instant  of  reception  tj.  If  each  observation  is  labeled  by  the  subscript 
of  the  time  at  the  instant  of  reception,  then  the  propagation  time  for  the  ith  observation  is  tpj,  where 

The  iteration  is  begun  with  tpj  = 0 on  the  right,  and  rs(tj)  is  calculated  by  the  interpolation  formula.  The 
result  is  a nonzero  value  for  tpj,  which  is  used  to  calculate  F^ftj  - tpj)  for  the  next  round.  Wlien  tpj  changes 
by  less  than  one  microsecond,  the  procedure  is  terminated.  The  range  for  the  time  tj  is  therefore 

Pi  = [(Xsfti  - tpi)-  Xg(ti))2  + (ys(ti  - tpi)  - yg(fi))^  + (zs(ti  - tpi)  - Zgfti))2]  (14) 

and  the  range  difference  calculated  for  the  observation  Oj  at  tj  is 


Pi  -Pi-i- 


One  subtle  point  that  must  be  mentioned  is  the  following.  The  satellite  position  is  calculated  on  the 
basis  of  an  earth  fixed  coordinate  system  (one  which  rotates  with  the  earth).  In  this  frame,  any  location 
which  is  fixed  to  the  earth  has  time  invariant  coordinates.  It  is  clear,  however,  that  the  signal  transmitted 
from  the  satellite  travels  in  inertial  space  not  in  a frame  rotating  with  the  earth.  Thus  each  time  Fg  is 
calculated,  a small  correction  must  be  made  to  the  station  position  to  account  for  the  station  motion  during 
the  propagation  interval.  Given  that  the  angular  velocity  of  the  earth  is  tOg  = 7.2921 15855  X 10"^  radians/ 
second,  in  a time  interval  tp  the  earth  will  have  rotated  through  an  angle  AX  = cOgtp.  With  the  station 
location  as  defined  by  Equations  (10),  the  corrections  to  Xg  and  yg  due  to  AX  (Zg  is  independent  of  X)  are 


AXg  = - (AX  sin  Xg)(0  + hg)  cos  0g 
Ayg  = (AX  cos  Xg)(0  + hg)  cos  </ig 

The  corrected  station  positions  are  therefore  Xg  = Xg  + Axg,  y^  = yg  + Ayg  and  Zg. 
ground  case  the  maximum  station  shift  is  about  7 meters. 


(15) 


In  the  NAVSAT-to- 


A small  correction  to  the  range  difference  (Equation  (15))  is  added  to  compensate  for  the  effect  of 
tropospheric  refraction.  The  Hopfield  Model'*  is  used  to  calculate  the  range  correction,  RH(t),  based  on  zenith 
angle  of  the  NAVSAT  as  seen  from  the  station,  and  local  temperature,  pressure,  and  humidity.  The  range 
correction  is  scaled  by  a factor  1 + Cr,  which  is  adjusted  in  the  solution  for  each  pass  independently.  The 
correction  for  each  observation  O,  is 

(1 +CR)[RH(ti)-RH(ti-i)]  (16) 


with  Cr  initially  set  to  zero. 


Also  included  in  the  solution  is  a frequency  bias  parameter  fn  - fb-  The  nature  of  this  term  is 
explained  in  the  previous  section;  its  form  is 


c(fn  - fb) 
f J + fo  + fb 


(ti-ti-l). 


(17) 


with  the  sum  fp  - fb  initially  set  to  zero. 


Summing  Equations  (IS),  (16),  and  (17)  provides  the  expression  for  Cj  used  to  form  the  observation 
matrix  D. 

. c(f„-fb) 

Ci  = (Pi-Pi.i)  + (l+CR)[RH(ti)-RH(ti_i)]  + , ^ (ti-ti_i)  (18) 

>1  + '0  ’b 

Differentiation  of  Equation  (18)  with  respect  to  the  parameters  Xg,  yg,  Zg,  Cr,  fb  forms  the  following  five 
partial  derivatives;  fp  is  assumed  zero. 


-^=  p''(ti_i)|X5(tj_,  - tpi.i)  -Xg(ti.i)]  -p-'(ti)(Xs(ti  - tpj)  - Xg(ti)] 

aCj 

= P''(ti.i )[ys(ti.i  - tpj.i) - yg(ti., )]  - p-'(ti)(ys(tj  - tpi)  - yg(ti)] 
aCj 


3Cj 

— =R„(ti)-R„(ti.,) 


8fb 


= -c(ti-tj_,) 


1^1  1^0 

(fi  + fo  + fb)^ 


-c(ti- Vl) 
1^0 


(19) 


When  these  are  evaluated  at  the  time  of  each  observation,  a matrbe  A is  formed  where  each  row  corresponds 
to  one  observation  time. 
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ac, 

ac, 

ac, 

ac. 

ac. 

dXgi 

9ygi 

3Cr, 

3'bi 

ac2 

ac2 

acj 

ac2 

ac2 

axg2 

3yg2 

azg2 

acR2 

9fb2 

ac„ 

ac„ 

aCn 

aCn  ac„ 

^ygn 

Corresponding  to  each  observation,  Oj,  there  is  a computed  range  difference,  Cj,  and  a linear  correc- 
tion consisting  of  partial  derivatives  from  a row  of  A and  incremental  variations  in  the  five  parameters 
previously  listed  and  henceforth  denoted  by  X. 

X = ISXj,  6yg  6Zp  6Cr  fiffe^. 


Therefore  in  matrbt  form 


O = C + AX, 

which  can  be  rearranged  by  defining  the  matrix  D: 

D = 0-C  = AX.  t:o) 

MATRIX  MANIPULATIONS 

Each  satellite  pass  observed  by  a ground  station  will  produce  a system  of  linear  equations  like  those 
represented  in  Equation  (20).  Solution  of  Equation  (20)  for  the  vector  X will  in  principle  provide  the  cor- 
rections needed  to  establish  the  true  values  for  Xg,  yg,  Zg,  Cr,  and  fj,.  The  matrix  manipulations  necessary 
to  obtain  the  least-squares  solution  arc  described  next. 

Each  observation  Oj  which  makes  up  Equation  (20)  is  assigned  a variance  which  is  dependent 
upon  the  duration  of  that  particular  Doppler  count.  The  expression  adopted  is  the  following: 

®obs  ~ over)"'  (seconds)"' 

The  reciprocal  of  each  individual  variance  is  assigned  to  the  diagonal  of  a square  weight  matrix  W.  The  off 
diagonal  elements  are  set  to  zero  while  the  ith  diagonal  element  is  Wjj  = 

The  first  operation  on  Equation  (20)  is  to  premultiply  by  A^W 

A^WAX  = A^WD 


13 


Figure  4.  Radial,  Along  Track,  and  Cross  Track  Unit  Vectors 


This  produces  the  normal  equations  denoted  by  the  square  matrix  B = A^WA.  The  matrix  product  on  the 
right-hand  side  is  E = ATWD.  The  above  expression  is  therefore  simplified  to 

BX  = E (21) 

If  a solution  for  X were  attempted  at  this  point,  the  result  would  be  indeterminate.  The  geometry 
involved  in  a single  pass  will  not  permit  a three-dimensional  solution  for  station  position.  However,  a two- 
dimensional  solution  is  possible  and  can  be  obtained  from  the  present  formulation  by  rotation  to  another 
coordinate  system.  Tlie  coordinate  frame  to  choose  is  one  which  will  isolate  the  geometric  ambiguity  to  a 
single  coordinate.  Such  a system  is  denoted  by  u,  v,  w (radial,  along  track,  cross  track),  illustrated  in 
Figure  4.  The  cross  track  direction,  w,  is  the  geometrically  weak  coordinate,  which  can  be  deweighted, 
forcing  the  least-squares  procedure  to  ignore  it  as  a possible  degree  of  freedom.  The  transf.>rmation  matrix 
R,  evaluated  at  the  time  of  minimum  slant  range,  operates  on  the  matrices  in  Equation  (21)  changing  the 
solution  vector  X to  Y,  where 


Y=  |6u6vSw5Cr  Sf^l'’’. 
14 


The  matrix  R is  derived  as  follows.  The  vector  to  the  NAVSAT  r^  is  as  given  in  Equation  (11),  while 
the  velocity  vector  is  r^ 

h = + VsY  + ZsZ. 

and 

The  three  unit  vectors  defining  the  new  coordinate  system  are 

u = Fj  r^  radial, 

w = (FjXFj)  ^ |rsXr,|  cross  track,  and 
V = wxu  along  track 

These  definitions  show  that  w is  perpendicular  to  the  plane  containing  the  radial  and  velocity  components 
of  the  satellite. 


The  rotation  matrix  R which  transforms  any  vector  X(x,  y.  z)  into  a vector  Y(u,  v,  w)  by  the  equation 
Y = RX  consists  of  the  following  elements: 


^21  = - ’‘s^-s)^s  - (’‘sVs  - ysXs)ysl 

R22  = [(’‘si's  - ysXs)Xs  - (Ys^-s  - ^sys)Zsl  - 
^23  = Kys^-s  “ ='-sys)ys  -(^s^s  - XsZs)Xsl  -i-  l(rsX^)xrsi 
R31  =(ysZs-Zsy,)H-|rsXr5l 
R32  = (^s^s  - ’‘s^-s)  IfsXrJ 

R33  = Kys-ys’‘s)^iwi. 

The  parameters  for  tropospheric  refraction  and  frequency  bias  are  not  altered  by  the  rotation,  thus  the  full 
rotation  matrix  is 


R = 


Rii 

Ri2 

Ri3 

0 

0 

R21 

R22 

R23 

0 

0 

R3, 

R32 

R33 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

1 

After  rotation,  a matrix  P containing  np  < 5 a-priori  variances  is  added  to  appropriately  initialize 
the  rotated  covariance  matrbe  by  deweighting  the  out-of-plane  coordinate  w.  The  rotation  is  performed  by 
pre-  and  postmultiplying  the  covariance  matrw  by  R and  R"’ , and  premultiplying  the  E matrix  by  R. 

RBR'lY  = RE 
15 
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Now  P can  be  added  to  the  rotated  covariance  matrix  resulting  in 

(RBR'l  + P)Y  = RE 

The  solution  vector  Y can  be  freed  by  premultiplying  both  sides  by  (RBR"'  + P)"'  leaving 

Y = (RBR-1  +P)-1  RE. 


In  terms  of  A and  W this  is 


Y = (RA^WAR-'  + P)-'  RA^WD 


(22) 


EDITING  PROCEDURE 

The  observations  belonging  to  each  pass  are  edited  independently  of  all  other  passes,  first  by  absolute 
tolerance  tests  and  then  by  evaluation  of  the  residuals  after  solution  of  Equation  (22).  The  adjusted  D 
matrix,  or  residual  matrix  is  then  calculated  to  establish  the  noise  level.  If  any  element  of  Dg^j  exceeds  the 
established  tolerance,  that  observation  is  eliminated  from  the  next  iteration  of  the  solution.  The  parameters 
which  are  used  for  the  editing  decisions  are  the  following: 

2 

"obs  “ observation  variance  which  appears  in  the  W matrix.  The  value  will  change  with  each 
edit  cycle. 

— The  initial  value  of  as  it  appears  on  the  observation  file. 

»?n,  - The  multiplication  factor  used  to  change  a^bs  from  one  edit  cycle  to  the  next.  Initially 
Vn,  = 1 and  consequently  = a,. 

T),  - A tolerance  factor  indicating  the  maximum  deviation  allowed  before  an  observation  is 
rejected.  Initially  rjj  » 1. 

Zg  - The  zenith  angle  tolerance.  Any  observations  obtained  at  zenith  angles  > are  rejected. 

- If  any  elements  of  the  D matrix  exceed  t?3,  they  are  rejected. 

One  edit  cycle  (which  is  concurrent  with  a solution  cycle)  begins  by  comparing  each  element  of  D 
with  t?3  and  each  zenith  angle  with  z^.  If  z > Zg  or  if  IDjl  > the  observation  is  eliminated  from  further 
consideration.  The  number  of  observations  remaining  in  a pass  after  these  tests  is  ng.  The  steps  which 
follow  are  numbered  for  convenience. 


I . If  IDjl  > Tj,T)p,a,,  the  observation  is  eliminated,  and  the  number  of  observations  lemaining  is 
decremented  by  one. 


2.  If  no  observations  are  rejected  in  Step  I of  this  cycle,  the  processing  of  this  pass  is  completed.  If 
this  is  the  first  cycle  through  the  editing  procedure,  continue  on  to  Step  3. 


3.  Recalculate  tj,  using  the  empirical  form 


n, 


0 


ng  + 50 


+ 2 


[—) 

\"r/ ' 


4.  Form  the  matrix  (Equation  (20)),  using  only  the  observations  remaining  and  continue  through  the 
manipulations  to  Equation  (22). 


5.  Calculate  the  parameter  V using  the  following  formula 


V=  Y — 

7 1 


then  use  V to  calculate  the  signal  to  noise 


^ V - E'^X 
N, 


Uq  + np  - 5 


Finally  recalculate  using  = (S/N)nn,  where  the  rj^  on  the  right  is  the  old  value,  and  t?m  on  the  left  is 
the  new  value. 

6.  Compute  the  adjusted  D matrix,  Dadj,  where 

D,dj  = D-AX 

Since  AX  is  the  linear  approximation  to  D (Equation  (20)),  the  difference  of  these  should  ideally  be  zero. 
However,  the  noise  in  Oj  will  prevent  a null  result.  What  is  desired  is  that  the  elements  of  the  matrix  be 
zero  mean,  white  Gaussian  noise.  The  variance  of  this  noise  is  the  quantity  of  interest  and  establishes  the 
performance  of  the  NAVPAC  hardware. 


7.  The  last  step  is  a return  to  Step  1 to  begin  the  next  edit  cycle.  In  all  cycles  other  than  the  first, 
the  matrix  is  used  in  place  of  the  D matrix  in  Steps  1 through  5. 


EDITING  RATIONALE 

The  editing  procedure  outlined  above  has  accreted  partly  empirically  and  partly  theoretically  over 
several  years  of  satellite  processing  at  NSWC/DL.  The  present  form  adaptively  adjusts  itself  to  the  peculiari- 
ties of  each  pass.  The  adjustments  are  applied  by  the  parameters  r?,  and  Tjp, , which  are  recalculated  in  each 
edit  cycle. 

The  product  (rj^o,)^  is  the  observation  variance  At  first  T?ni  . which  makes  the  observation 
variance  just  the  initial  estimate  a^.  Each  succeeding  edit  cycle  should  see  decrease,  thus  restricting  the 
number  of  observations  which  pass  Step  I . The  parameter »?,  multiplies  by  some  factor  greater  than 
two,  and  thus  allows  observations  to  pass  Step  I if  they  deviate  less  than  from  the  mean.  The 
calculation  of  ry,  in  Step  3 depends  upon  the  number  of  observations,  nQ,  in  the  pass,  and  the  number 
remaining  at  the  current  edit  cycle  step  n, . A NAVSAT  pass  as  seen  from  the  ground  will  rarely  contain  30 
observations;  thus  when  n,  = n^  the  value  of  ry,  < 2.375.  As  observations  are  eliminated,  ty,  will  increase 
slightly. 

The  calculation  of  ry^  measures  how  well  the  present  solution  has  performed  as  compared  with  the 
previous.  This  is  expressed  by  the  signal  to  noise  in  Step  5.  Contrary  to  most  definitions  of  signal  to  noise, 
the  object  here  is  to  make  the  ratio  unity.  The  signal  in  this  case  represents  dynamics  which  have  not  been 
removed  by  the  least-squares  model.  If  the  model  truly  represents  the  observations,  then  only  white 
Gaussian  noise  can  be  left. 


Following  OToole* , the  value  for  the  multiplying  factor  tj^  is  adjusted  so  that  the  signal  to  noise 
squared,  defined  by 


"L  J’ 

is  approximately  unity.  In  the  above,  n is  the  sum  Oq  + Up  - 5,  and  C'  is  defined  below. 

To  show  that  this  form  of  (S/N)^  is  the  same  as  that  given  in  the  preceding  section,  it  is  necessary  to 
reintroduce  the  computed  vector  C,  the  observation  vector  O,  plus  the  weight  matrix  W.  A new  computed 
vector  C',  containing  the  computed  values  after  the  solution,  is  defined  by 

C'  = C + AX 

An  expanded  form  of  equation  (21)  is  also  needed 

BX  = E = a’*'w(0  - C),  (23) 


along  with  the  variance  V rewritten  in  matrix  form 


V = ! (0  - C)'^(0  - C). 


It  is  assumed  that  is  the  same  for  all  observations. 

In  order  to  procede,  the  quantity  (0  - C')^(0  - C')  must  be  expanded  by  substituting  for  C'. 


(O  - C')'’'(0  - C')  = (O  - C)^(0  - C)  - (0  - C)^ AX  - (AX)^(0  - C)  + (AX)”^ AX 


Dividing  the  above  by  allows  the  weight  matrix  to  be  introduced.  For  simplicity  it  would  be  de- 

fined as  W = (l/(T?n,a,)^)l. 


(O-C')^(O-C') 


(O  - C)'^(0  - C) 


- 2(AX)'^W(0  - C)  + X'^a'^WAX 


Substituting  for  the  variance  and  using  (23)  gives 


= V-2X'^E  + X^E, 


which  can  be  simplified  to  the  desired  form: 


- (V-E^X) 
n 
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The  signal  to  noise  is  used  to  scale  the  multiplying  parameter  TJm,  which  operates  on  o,  to  set  a„bs 
Step  1 of  the  editing  cycle.  After  several  cycles,  this  feedback  allows  the  number  of  observations  rejected 
to  reach  a steady  state  and  makes  the  tolerance  testin'-  less  sensitive  to  the  value  chosen  for  o,. 


EXAMPLES  OF  ADJUSTED  RESIDUALS 

Examples  of  adjusted  residuals  are  shown  in  Figures  5 through  16.  Adjusted  D matrices  obtained  by 
processing  real  data  are  shown  in  Figures  5,  8,  and  1 1.  For  comparison.  Figure  14  is  computer-generated 
zero  mean  Gaussian  white  noise  with  unit  standard  deviation.  Figure  1 5 illustrates  a histogram  generated 
from  this  sample  of  white  noise  with  the  mean  and  standard  deviation  calculated  from  the  histogram.  The 
dashed  line  is  a Cuassian  function  having  the  same  mean,  standard  deviation,  and  maximum  as  the  histogram. 
Figure  16  is  the  autocorrelation  function  calculated  from  this  same  sample.  If  the  computer-generated  noise 
were  perfectly  white,  there  would  be  no  correlation  and  the  autocorrelation  function  would  consist  of  an 
impulse  at  zero  lag.  Figures  5 through  13  are  plots  similar  to  those  described  above  but  obtained  from  the 
single  pass  NAVPAC  solution  and  editing  procedure  outlined  here. 

An  interesting  feature  of  the  autocorrelation  functions  is  that  they  are  distinctly  different.  Figure  7 
shows  two  definite  secondary  correlations  at  lag  steps  9 and  16.  Figure  10  docs  not  distinguish  itself  with 
obvious  correlations  but  is  more  like  the  white  noise  result  than  the  other  two.  The  function  plotted  as 
Figure  13  is  intermediate  between  the  first  two  showing  possible  regions  of  correlation  both  positive  and 
negative. 

The  different  autocorrelation  results  illustrated  seem  to  indicate  time  variable  phenomena  influencing 
the  data  acquisition  portion  of  NAVPAC.  Since  all  three  cases  were  processed  in  the  same  way  using  the 
same  model,  any  variation  in  the  results  ist  be  due  to  either  propagation  effects  or  peculiar  time  variable 
correlations  introduced  by  the  NAVPAC  electronics.  These  autocorrelation  functions  are  a quick  visual 
guide  to  the  quality  of  the  individual  pass.  Strong  correlations  indicate  a poor  fit  of  the  observed  data  to 
the  model;  weak  or  rapid  variations  indicate  a good  fit. 
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